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ABSTRACT 



The structure of the inner region of an advection-dominated accretion disk around 
O , 

a nonrotating black hole is explored by applying asymptotic analysis in the region 
just outside the event horizon. We assume that the viscous transport is described by 
the standard Shakura-Sunyaev prescription throughout the disk, including the inner 
region close to the horizon. One of our goals is to explore the self-consistency of this 
assumption by analyzing the causality of the viscous transport near the black hole. The 
effects of general relativity are incorporated in an approximate manner by utilizing a 
pseudo-Newtonian gravitational potential. Analysis of the conservation equations yields 
unique asymptotic forms for the behaviors of the radial inflow velocity, the density, 
the sound speed, and the angular velocity. The specific behaviors are determined by 
three quantities; namely, the accreted specific energy, the accreted specific angular 
momentum, and the accreted specific entropy. The additional requirement of passage 
through a sonic point further constrains the problem, leaving only two free parameters. 
Our detailed results confirm that the Shakura-Sunyaev viscosity yields a well-behaved 
flow structure in the inner region that satisfies the causality constraint. We also show 
that the velocity distribution predicted by our pseudo-Newtonian model agrees with 
general relativity in the vicinity of the horizon. The asymptotic expressions we derive 
therefore yield useful physical insight into the structure of advection-dominated disks, 
and they also provide convenient boundary conditions for the development of global 
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models via numerical integration of the conservation equations. Although we focus here 
on advection-dominated flows, the results we obtain are also applicable to disks that 
lose matter and energy, provided the loss rates become negligible close to the event 
horizon. 

Subject headings: accretion disks — hydrodynamics — black holes — methods: analyt- 
ical — general relativity 

1. INTRODUCTION 

The advection-dominated accretion flow (ADAF) model has received a great deal of attention 
as a possible explanation for the dynamics occurring in X-ray under luminous, radio-loud Active 
Galactic Nuclei (AGNs), which are thought to contain hot accretion disks (e.g., Narayan & Yi 
1994, 1995). Briefly, these models describe the dynamics of gas fed onto a black hole at very low 
(significantly sub-Eddington) accretion rates. Gas accreting at such a low rate is quite tenuous, 
and therefore the ion-electron Coulomb coupling timescale can exceed the timescale for accretion. 
Since the ions absorb most of the energy dissipated via viscosity, and the Coulomb coupling with 
the electrons is weak, the ions achieve a nearly virial temperature (Tj ~ 10 12 K) that greatly 
exceeds the electron temperature (T e ~ 10 9 K), unless plasma instabilities directly heat the electrons 
(Bisnovatyi-Kogan & Lovelace 1997). The tenuous nature of the gas therefore severely limits the 
radiative efficiency of the plasma, and consequently most of the thermal energy dissipated by 
viscosity is advected into the black hole in the form of hot protons, although outflows may carry 
some of this energy away. The resulting X-ray luminosity is far below the Eddington limit. 

In X-ray bright AGNs, the disk is thin and cool, and the accretion proceeds in a radiatively 
efficient manner (Narayan 2002). This explains the origin of the "big blue bump" in the spectra 
of typical Seyfert galaxies, along with the common occurrence of broad emission lines, appar- 
ently formed in the inner region between the last stable orbit and the event horizon. The high 
temperatures in ADAF disks preclude the formation of either the blue bump or the broad lines. 
From a theoretical point of view, the reason a given object chooses one mode of accretion over 
the other is not entirely clear. As Narayan (2002) points out, it may be possible to learn a great 
deal about how the flow changes character as a function of luminosity by studying transition ob- 
jects such as low-ionization nuclear emission regions (LINERs) and low-luminosity active galactic 
nuclei (LLAGNs). Ptak et al. (1998) suggest that the larger characteristic variability timescales 
observed in these sources compared with brighter Seyfert galaxies may indicate the presence of 
central ADAFs, which increase the size of the emission region relative to more efficient, thinner 
disks. This idea is supported by observations of the bright Seyfert 1 galaxy IC 4329A performed 
by Done, Madejski, & Zycki (2000) using ASCA and RXTE. They note that the iron line in this 
source is not as broad as that detected in more extreme cases such as MCG-6-30-15, where the cool 
disk apparently extends down to the last stable orbit. Based upon this observation, they conjecture 
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that in IC 4329 A the cool disk transitions into a central ADAF outside the last stable orbit, and 
consequently the inner region is too hot to produce the line emission. The physical nature of the 
accretion disk may also vary as a function of time in individual sources. For example, based on 
intensive RXTE observations of the galactic black hole candidate J1550-564, Wilson & Done (2001) 
propose that transitions between the low/hard and high/soft spectral states reflect the appearance 
and disappearance of an ADAF in the inner region, perhaps as a consequence of changes in the 
accretion rate. There are indications of similar behavior in some AGNs (e.g., Lu & Yu 1999). 

These observations make it clear that in order to unravel the complex global structure of the 
accretion disk, a complete understanding of the physical properties of ADAFs close to the event 
horizon is essential. One of the most important unresolved questions concerns the behavior of the 
torque in the inner region, where the material begins to plunge into the black hole. In the standard 
thin-disk model, it is usually assumed that the stress vanishes at the marginally stable orbit (e.g., 
Frank, King, & Raine 1985). Others have used causality arguments to suggest that the stress 
vanishes at the transonic point (e.g., Popham & Narayan 1992). However, several authors have 
argued based on the results of magnetohydrodynamical simulations that magnetic stresses are able 
to remove angular momentum from the material in the plunging region (e.g., Reynolds & Armitage 
2001; Hawley & Krolik 2001; Agol & Krolik 2000; Gammie 1999). The only model-independent 
statement that can be made with absolute certainty is that the horizon itself cannot support a 
shear stress. Hence we view this as the most conservative possible hypothesis. The proposition 
that the viscous torque actually vanishes at some radius outside the horizon is inevitably model- 
dependent. Following Narayan, Kato, Sz Honma (1997), Yuan (2001), and Yuan et al. (2000), we 
shall therefore assume here that the torque vanishes at the horizon, and ask whether self-consistent 
ADAF models with this property can be constructed based upon the standard Shakura-Sunyaev 
viscosity prescription. Our basic approach is to explore the associated disk structure using rigorous 
asymptotic analysis. The validity of the dynamical results is evaluated by examining the causality 
of the viscous transport near the horizon, where all signals must propagate into the black hole. We 
also compare the velocity distribution in the vicinity of the horizon with the predictions of general 
relativity. Based upon these considerations, we argue that a self-consistent Shakura-Sunyaev flow 
can exist in the inner region, and we present detailed asymptotic solutions for the disk structure. 

The remainder of the paper is organized as follows. In § 2 we briefly review the sequence 
of models developed to describe ADAF disks. The conservation equations for one-dimensional 
ADAF disks are discussed in § 3. The appropriate inner boundary conditions for ADAF disks 
are derived in § 4 by employing asymptotic analysis of the conservation equations. In § 5 the 
boundary conditions are used to obtain global solutions for the physical quantities by numerically 
integrating the conservation equations. The exact numerical solutions obtained are compared with 
the asymptotic expressions and also with solutions previously presented in the literature. The 
implications of our results for the structure of advection-dominated accretion disks around black 
holes are discussed in S 6. 
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2. DISK MODELS 

The structure of advection-dominated disks has been explored using a variety of theoretical and 
computational approaches. The initial one-dimensional, self-similar models developed by Narayan 
& Yi (1994, 1995) incorporated Newtonian gravity. In later works, the assumption of self-similarity 
was relaxed and new one-dimensional models were developed based on a complete set of conservation 
equations (Narayan, Kato, & Honma 1997; Chen, Abramowicz, & Lasota 1997). These models 
utilized a pseudo-Newtonian gravitational potential in order to approximate the effects of general 
relativity (Paczyhski & Wiita 1980; Abramowicz, Calvani, & Nobili 1980). Conservation of mass, 
energy, and angular momentum is expressed using a set of coupled differential equations that are 
integrated to obtain the inflow velocity v, the sound speed a, and the angular velocity as functions 
of radius. The solution of these equations requires a sufficient number of boundary conditions, 
imposed either at the sonic point, the black hole horizon, or the outer radius of the computational 
domain. The sonic point is a critical point for the flow, and therefore the integration must be 
divided into two regions. One typically solves the equations by starting at the sonic point with 
values that satisfy the critical conditions, and then integrating the differential equations inward 
towards the horizon, and outward toward larger radii, where the disk is eventually expected to 
become thin and cool. 

Because of the singularity at the event horizon introduced by the pseudo-Newtonian potential, 
direct integration of the conservations equations towards the horizon using an explicit method 
such as a Runga-Kutta algorithm cannot be used to obtain smooth global solutions. On the other 
hand, explicit integration away from the horizon towards the critical point is stable, and therefore 
it represents a more convenient approach to the problem, if the boundary conditions close to the 
black hole event horizon can be specified. Another alternative is to employ a relaxation algorithm 
that is based on iteration of the numerical solution, with the goal of minimizing a global error 
parameter (Press et al. 1986). Each of these methods requires the availability of a suitable set of 
inner boundary conditions describing the physics of the gas close to the event horizon. However, 
the boundary conditions appropriate for this problem have not been presented previously in the 
literature. Motivated by the lack of this crucial information, in this paper we derive the exact 
asymptotic forms for the variation of the physical quantities close to the event horizon by employing 
asymptotic analysis based on the differential conservation equations. We shall specialize to the case 
of one-dimensional, advection-dominated flow in the pseudo-Newtonian potential. The asymptotic 
results will be used to develop inner boundary conditions that facilitate the integration of the 
conservation equations. 

3. CONSERVATION EQUATIONS 

We shall focus on the structure of steady, advection-dominated accretion disks. These disks 
accrete at well below the Eddington rate, and are therefore so tenuous that radiative cooling is 
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inefficient. Hence advection-dominated disks are nearly virially hot, and are essentially collisionless. 
Moreover, the electrons and the protons are likely to possess distinct temperatures due to the low 
rate of Coulomb interactions between these species. In this situation, the protons absorb most of 
the energy dissipated via viscosity, and consequently they possess a much higher energy density 
than the electrons. Height-integrated structure equations for such disks have been derived by 
Abramowicz et al. (1988). 

In our approach to modeling the disk structure, we will incorporate the effects of general 
relativity in an approximate manner by expressing the gravitational potential per unit mass using 
the pseudo-Newtonian form (Paczyhski & Wiita 1980) 



where r s = 2GM/c 2 is the Schwarzschild radius for a black hole of mass M. This potential gives 
correct results for the location of the horizon, the radius of marginal stability, and the radius of the 
marginally bound orbit around a nonrotating black hole (Paczyhski & Wiita 1980; Abramowicz, 
Calvani, & Nobili 1980). Utilization of equation (1) is convenient because it facilitates a semi- 
classical approach to the problem that simplifies the analysis considerably, while maintaining good 
agreement with fully relativistic calculations. For this reason, the pseudo-Newtonian potential has 
been adopted by a number of authors in their investigations of accretion onto Schwarzschild black 
holes (e.g., Matsumoto et al. 1984; Abramowicz et al. 1988; Chen et al. 1997; Narayan et al. 1997; 
Hawley & Krolik 2001, 2002; Yuan 1999; Yuan et al. 2000; Reynolds & Armitage 2001). We show 
in Appendix A that the relativistically correct energy equation for a particle freely falling from rest 
at infinity in the Schwarzschild metric can be written as 



where v r and v v denote the radial and azimuthal components of the four-velocity, respectively, 
and 3>(r) is given by equation (1). The Newtonian appearance of equation (2) can be viewed as 
one of the primary motivations for introducing the pseudo-Newtonian potential ^(r), although one 
needs to keep in mind that the dynamical quantities v r and v v introduced in equation (2) are four- 
velocities rather than conventional velocities. We will return to this point later in our discussion of 
the dynamical results. 



In our subsequent analysis, we will adopt the "perfect ADAF" approximation, meaning that 
(for now) we shall completely neglect the escape of energy and matter from the disk. Modifications 
associated with the relaxation of this assumption will be discussed in § 6.4. In the one-dimensional, 
steady state ADAF scenario, three quantities are conserved in the flow, namely, the accretion rate 




(2) 



3.1. Transport Rates 



M = 2vrr • 2H ■ pv , 



(3) 
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the angular momentum transport rate 

j = Mr 2 Q-g, (4) 

and the energy transport rate 

E = -gn + M^w 2 + ^v 2 + ^-^ + ^j , (5) 

where p is the mass density, v is the radial velocity (defined to be positive for inflow), Q is the 
angular velocity, w = r f2 is the azimuthal velocity, U is the internal energy density, P is the 
gas pressure, and g is the torque. All quantities represent vertical averages over the disk half- 
thickness H. We shall assume that the ratio of specifics heats, 7 = (U + P)/U, maintains a 
constant value throughout the flow. Note that the transport rates M, J, and E are all defined to 
be positive for inflow. Since these quantities are conserved, they represent the rates at which mass, 
angular momentum, and energy are swallowed by the black hole. Athough we refer to v and w as 
"velocities," we shall see later that based on their asymptotic behavior close to the event horizon, 
these quantities are actually more correctly interpreted as four-velocities. 



3.2. Momentum Equations 

The expressions for the transport rates M, J, and E are supplemented by equations describ- 
ing the conservation of the three components of momentum. In most of the one-dimensional disk 
models, the thickness of the disk is computed using the assumption of vertical hydrostatic equilib- 
rium. While this assumption is probably well-satisfied in the regions of the disk that have subsonic 
radial velocities, it may not be very accurate near the horizon because the flow is supersonic and 
practically in free-fall there. Despite this, the hydrostatic assumption is routinely used to describe 
the entire disk, including the supersonic region (e.g., Chen et al. 1997; Narayan et al. 1997; Chen 
& Taam 1993). Since one of our motivations is to develop a consistent set of inner boundary condi- 
tions applicable to "standard" ADAF models, we shall assume for now that the disk half-thickness 
is given by the usual hydrostatic prescription, 

H(r) = b £ , (6) 

where 60 is a dimensionless constant of order unity that depends on the details of the vertical 
averaging (Abramowicz et al. 1988), 

fp\ 1/2 

« W S (-) (7) 

represents the isothermal sound speed, and QkW denotes the Keplerian angular velocity of matter 
in a circular orbit at radius r in the pseudo-Newtonian potential (eq. [1]), defined by 

n2M GM 1# , , 

r{r — r s Y r ar 
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In § 6.2, we will discuss how our results would be modified if the assumption of vertical hydrostatic 
equilibrium in the supersonic region were replaced with radial free-fall. 

In a steady state, the comoving radial acceleration rate in the frame of the accreting gas is 
expressed by 

Dv - dv - 1 dP d ® r n 2 

Dt dr p dr dr 
and angular momentum transport is treated by relating the torque Q to the gradient of the angular 
velocity O using the fundamental formula (e.g., Frank, King, & Raine 1985) 

Q = -AnrH pur 2 , (10) 
dr 

where v is the kinematic viscosity. 



3.3. Viscosity and Torque 

We shall adopt the Shakura-Sunyaev (1973) prescription for the kinematic viscosity, 

a a 2 

where a is a constant. While this is certainly a reasonable approach in the subsonic region of 
the flow, the validity of this prescription in the supersonic region depends on the microphysical 
mechanism responsible for generating the torque in a given situation. In particular, if the angular 
momentum is transferred via turbulent blobs of fluid, and the turbulence is subsonic, then the 
associated torque should vanish in the supersonic region as a consequence of causality considerations 
(e.g., Kato 1994; Narayan 1992; Popham & Narayan 1992). However, it is not clear a priori whether 
fluid turbulence, particles, or magnetic fields (or some combination of these effects) transports the 
angular momentum in ADAF disks. If the transport occurs via particles or fields, then the causality 
argument mentioned above does not apply, and torques can be generated even in the supersonic 
region of the flow between the sonic point and the horizon, although the torque must certainly 
vanish at the horizon itself in keeping with general relativistic considerations. In fact, simulations 
performed by Reynolds & Armitage (2001) suggest that magnetic fields are able to transfer angular 
momentum from material in the plunging region to material in the outer disk. Following Narayan et 
al. 1997, we will therefore assume here that the angular momentum is transferred by some generic 
particle/magnetohydrodynamical mechanism, characterized by an effective value for the adiabatic 
index 7 in the range 4/3 < 7 < 5/3. In this case, we can safely adopt the standard a-prescription 
for the kinematic viscosity v given by equation (11). The self-consistency of this approach will be 
evaluated in § 6.3, where the causality of the flow near the horizon is examined. 
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3.4. Entropy and Internal Energy 

Since we are neglecting the escape of energy from the disk, the comoving rate of change of the 
internal energy density U can be written in the frame of the gas as 

DU dU U dp ■ 

-prr = -v— = -7 — v— + (7 v iscous , (12) 
Dt dr p dr 

where 

is the viscous energy dissipation rate per unit volume. The negative signs appearing in equation (12) 
reflect the fact that v has been defined to be positive for inflow. Combining equations (11), (12), 
and (13), we can rewrite the internal energy equation as the entropy equation 

d, fU\ U viscous a( 7 -l)r 2 f <Kl\ 2 

This equation demonstrates that the flow approaches a purely adiabatic behavior (U oc p 1 ) wherever 
the viscous dissipation rate U w - lscons /U vanishes. If the gas is in local thermodynamic equilibrium, 
then the viscous heating is a quasi-static process, and in this case the flow is isentropic wherever 
the dissipation vanishes. 

We shall find it convenient to express the variation of the isothermal sound speed, a, using the 
"entropy function," 

2+1 

Kir) . ^ . (15) 

To understand the physical significance of K, we can combine equations (3), (6), (7), and (15) to 
show that 

K 1-1 oc — , (16) 

which establishes that K remains constant in regions of the flow unaffected by dissipation. In 
particular, if the gas is in local thermodynamic equilibrium, then we can use equation (16) to 
demonstrate that the value of K is related to the entropy per particle S by (Reif 1965) 

S = klnK + c , (17) 

where cq is a constant that depends only the composition of the gas, but is independent of its state. 
Note that the relation between K and S in equation (17) may be violated in an ADAF because 
the gas is collisionless, and it is uncertain whether collective processes can establish a Maxwell- 
Boltzmann distribution. However, in any case K itself is unambiguously defined by equation (15). 
By comparing equations (14) and (16), we can show that the radial derivative of K is given by 

dlriK ^visCOUS _ < ^fi( < ^}\ 2 f-|Q\ 

V dr ~ (7- 1)U ~ ft K \dr J ' { j 
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This result confirms that K remains constant in regions of the flow that are not subject to dis- 
sipation. We will employ equation (18) in § 4.4, where we derive the asymptotic solution for the 
variation of the entropy function close to the horizon. 



4. ASYMPTOTIC ANALYSIS 

The conservation equations presented in § 3 can be solved as a coupled set to determine 
the radial profiles of the physical quantities v, S7, a, p, P, and H. Due to the divergence of 
the pseudo-Newtonian potential as r — > r s , the event horizon is a regular singular point of the 
conservation equations governing the disk structure. It is therefore possible to develop Frobenius- 
style expansions of the physical quantities around the point r = r s . Rather than developing 
complete series solutions for the variables, we will employ asymptotic analysis to determine the 
dominant behaviors as r — > r s . As we demonstrate below, this information can be used to derive 
explicit boundary conditions applicable very close to the horizon. 



4.1. Stress Boundary Condition 

One of the fundamental boundary conditions for black hole accretion is that the viscous shear 
stress S must vanish as r — > r s , because particles at different radii become causally disconnected 
from each other. The shear stress (force per unit area) is related to the torque Q by 

E = -""■* -SPff- <19> 

and consequently the vanishing of the stress on the horizon implies that Q = there as well. We 
therefore conclude based on equation (4) that 

hm n(r) = n = -f-, (20) 

which we can rewrite in terms of the azimuthal velocity w = r as 

lim w = — , (21) 

where 

f »4 (22) 

is the specific angular momentum of the material entering the black hole. Hence w approaches a 
finite value at the event horizon. Let us consider the physical implications of this result. Recall 
that the azimuthal velocity of a freely-falling particle as measured by a stationary observer 
outside a Schwarzschild black hole vanishes at the horizon (Shapiro & Teukolsky 1983). Hence, w 



-10- 



cannot represent the true azimuthal velocity measured by a static observer in the region close to 
the horizon. In Appendix A, we demonstrate that w is actually equal to the azimuthal component 
of the particle's four-velocity, v v , and consequently w possesses a finite value at the horizon. 



4.2. Radial Velocity 

Our result for the asymptotic variation of f2 close to the horizon (eq. [20]) can be used to 
derive the corresponding behavior of the radial inflow velocity v as r — > r s . By using equation (4) 
to eliminate the torque in equation (5), we can rewrite the energy transport equation as 

E = jn + M(-v 2 --r 2 n 2 + ^—a 2 --^-) , (23) 
\ 2 2 7 — 1 r — r s / 

where we have also substituted for the potential using equation 1. The flow into the black hole 
must be supersonic at the horizon since the radial velocity approaches the speed of light there. 
Hence !)>«asr-»r s . Since we have determined that the angular velocity O approaches a finite 
value at the horizon, we can conclude based on equation (23) that the radial velocity approaches 
the free-fall velocity vg(r), i.e., 

,. v ( r ) i , x (2GM\ l/2 

lim -y- = 1 , veir) = . 24 

• Vff(r) \r-r s J 



i — >r, 



This implies that close to the horizon, v formally exceeds c, and therefore it cannot represent 
the actual radial velocity v r measured by a static local observer in the Schwarzschild metric. In 
Appendix A we demonstrate that close to the horizon, v is actually equal to the radial component 
of the four- velocity, v r , for a freely falling particle. 



4.3. Angular Momentum and Torque 

The vanishing of the stress at the horizon ensures that the disk experiences differential rotation 
with dd/dr < at all radii. We can use this observation along with the requirement that — ^ S7o 
as r — > r s to develop the leading-order behavior of the Frobenius expansion for fi(r) about r = r s . 
Since r = r s is a regular singular point of the conservation equations, we can in general write the 
asymptotic behavior of Q(r) close to the horizon as 

n(r) = n -A(r-r s y , (25) 

where A and q are positive constants and we will use the symbol "=" to denote asymptotic equality 
at the horizon. Equation (25) is the simplest form that satisfies the requirements that — > Oo 
and dQ/dr < as r — > r s . The right-hand side of equation (25) represents the first two terms 
of the Frobenius expansion for f2(r), and q is the exponent of the solution about r = r s (Boyce 
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fe DiPrima 1977). We can constrain the value of q by examining the associated variation of the 
specific angular momentum, £ = r 2 Q,. Differentiation of £ with respect to radius yields 

^=2rfi + r 2 ?. (26) 
dr dr 

Substituting for using equation (25), we obtain the asymptotic relation 

AO 

^- = 2rn -2Ar(r-r s y-Aqr 2 (r-r s ) q - 1 . (27) 

Now, according to equation (4), £ = ( J + Q)/M, and therefore d£/dr > at the horizon since the 
torque Q vanishes there. This in turn implies that d£/dr has a finite, positive value at the horizon. 
Based on this constraint, we conclude that q > 1, since otherwise d£/dr would diverge to negative 
infinity at the horizon. 

The conclusion that q > 1 implies that dVi/dr approaches a. finite value as r — > r s . Furthermore, 
since v and fix each diverge as r — > r s , we can demonstrate based on equation (14) that the flow 
displays a purely adiabatic behavior close to the horizon (i.e., U oc p" 1 ). This explicitly confirms 
that the dissipation vanishes at the horizon, which is of course intuitively obvious since the stress 
vanishes there. According to equation (18), the entropy function K consequently approaches a 
finite value at the horizon, i.e., 

lim K{r) = K , (28) 

r^r s 

where k In Kq+cq represents the specific entropy of the particles entering the black hole (see eq [17]). 

We can build on our previous conclusions to further explore the asymptotic behavior of the 
specific angular momentum close to the horizon. Combining equations (3), (4), (10), (11), and 
(22), we can express the radial derivative of Q as 

dn _ v n K (£-£ ) 

i 2 2' V / 

dr a r L a z 

Using this result to substitute for d£ljdr in equation (26) yields a differential equation for £, 

d£ 2£ vn K (£-£ ) 



dr r a a 



2 



(30) 



As was pointed out earlier, £ = ( J + Q)/M, and therefore d£/dr > at the horizon since Q — > as 
r — > r s . It follows that the local behavior of £ close to the horizon must be of the general form 

l{r) =£ + B(r-r s f , (31) 

where B and f3 are positive constants. This represents the leading behavior of the Frobenius 
expansion for £(r) about r = r s , and (3 is the exponent of the solution. Paczyhski & Wiita (1980) 
imposed equation (31) as an ad hoc expression for the global variation of the specific angular 
momentum £(r), whereas we employ it only in the asymptotic limit, where its validity is fully 



- 12 - 



supported by the conservation equations. In fact, we shall see that this is not a very accurate 
expression for £ far from the event horizon in our application. Using equation (31) to substitute 
for £ on the right-hand side of equation (30) yields in the limit r — > r s 

Mm V -^B(r-r s r= 2 -^-l' Q , (32) 

r-*r s ar r s 

where 

e = Mm . (33) 

r^r g dr 

The constants B, (3, and £' are determined as follows. To evaluate the limit on the left-hand 
side of equation (32), we substitute for a using equation (15) and set K = K at the black hole 
horizon. We have already determined that the gas approaches radial free-fall as r — > r s , i.e., 
v — > [2GM/(r — Tg)] 1 / 2 . Using this information, we obtain 



2 1 / 2 GM 

hm 



2r 3 ( r 



r ^ r s aT 1 / 2 (r — r s ) 3 / 2 



7-1 



7+1 5(r-r s )^ = ^-^. (34) 

T Q 



In order to obtain a constant value on the left-hand side in the limit r — > r s , we must require that 
the exponents of (r — r s ) add to zero. This yields the result 

" = ^)' (35) 

In Table 1 we list the values of j3 obtained for several values of 7. Note that for 7 is the range 
4/3 < 7 < 5/3, we find that 1.36 > (3 > 1.25. Hence (3 exceeds unity for any physically acceptable 
equation of state. Using equation (31) to evaluate di/dr in the limit r — > r s therefore yields 

£' = lim 5 (3(r - ^ f' 1 = . (36) 

Hence we have proven that dl/dr = at the black hole event horizon. By balancing the values on 
the two sides of equation (34) in the limit r — > r s , it is straightforward to show that the constant 
B is given by 

Combining equations (31), (35), and (37), we can express the asymptotic solution for £ near the 
horizon as 

This relation gives the dominant asymptotic behavior of the specific angular momentum as a 
function of r for arbitrary values of a, 7, £q, and Kq. It is interesting to note that although 
Paczyhski & Wiita (1980) arbitrarily imposed equation (31) as a global expression for £(r), our 
numerical results for (3 are relatively close to the values they obtain. 
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Our asymptotic result for £{r) has two important implications. First, since we have found 
that dt/dr = at the event horizon, it follows from consideration of equation (27) that q = 1 and 
A = 2Qo/r s . Referring to equation (25), we conclude that the asymptotic solution for f2(r) is 
therefore given by 

n(r)=n -^(r-r s ) (39) 



in the vicinity of the horizon. Note in particular that dd/dr = — 2i7o/ r s a ^ the horizon, in 
contradiction to Narayan et al. (1997), who erroneously stated that d£l/dr = there. Second, 
since the torque Q is linearly related to £ via Q = Ml — J = M{£ — £o), we find that the asymptotic 
variation of the torque is given by 

nt , . a£ M ( 2 \ 1/2 ( Kl \^ . ,t±6. 



GM \rj \2r3 

Based on this expression, we conclude that the radial derivative of the torque vanishes at the horizon, 
i.e., 

lim % = . (41) 

r^r g dr 

This completely new boundary condition is one of the main results of the paper. The vanishing of 
the derivative dQ/dr at the horizon supplements the well-known boundary condition Q = 0. The 
physical interpretation of this new boundary condition and its effect on the structure of the global 
flow solutions will be discussed in 5 and 6. 



4.4. Entropy 

Our insights regarding the asymptotic behaviors off, £1, a, and v can be combined to ascertain 
the asymptotic variation of the entropy function K{r) close to the horizon. Substituting into 
equation (18) for dXuKjdr using the asymptotic free-fall velocity (eq. [24]) along with equation (8) 
for gives the leading behavior 

near to the horizon, where we have also used the fact that dQ/dr = —2 £lo/r s at r = r s . Integration 
of equation (42) with respect to radius yields for the asymptotic behavior of K(r) the solution 

(43) 

in the vicinity of the horizon. Note that K — > Kq at the horizon as required, although the radial 
dependence is very weak, reflecting the gradual disappearance of viscous dissipation as r — > r s . 
The flow is therefore essentially isentropic close to the horizon. 



K(r) = K 



5GM V 2 



1/2 



(r-r s f/* 
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4.5. Sound Speed and Inflow Velocity 

We have shown that v approaches the radial free-fall velocity vq = [2GM/(r — rg)] 1 / 2 near 
the horizon. While this result is certainly valid in the limit r ^ r s , we may strive to obtain a 
more precise formula for v by employing the insights we have obtained in our study of the local 
variations of the specific angular momentum £{r) = r 2 fi(r) and the entropy function K(r) close to 
the horizon. Based on equation (23), we can express the energy per unit mass transported through 
the disk as 



E v 2 I 2 £ £ 7 2 GM 

e = ^ = — -—^ + \ + —!— a 2 , 44 

M 2 2r 2 r 2 7-I 



where £q = J/M is the accreted specific angular momentum. Since E and M are conserved, it 
follows that eo represents the energy per unit mass swallowed by the black hole. 

Close to the horizon, the viscous dissipation vanishes, and K — ► Kq. We can therefore use 
equation (15) to express the variation of the isothermal sound speed a in the vicinity of the horizon 

as 

•W=H S (45) 
\ rv J 

We also know that l^^asr^r s . Using this condition along with equation (45), we can rewrite 
equation (44) as the asymptotic expression 

. v 2 £l 7 / r 2 v 2 \^ GM 
e ° = T + 2^ + 7^Up4J ~~s' (46) 
This nonlinear relation governs the variation of the inflow velocity v(r) close to the horizon. In 
general, it must be solved numerically to determine v for given values of £0, Kq, eo, and r. However, 
as an alternative to numerical root-finding, we can seek an asymptotic, analytical solution for v(r) 
by expanding equation (46) in the small parameter g(r), defined by 

v 2 (r) = v 2 (r)[l+g(r)} . (47) 

Very close to the horizon, the velocity approaches free-fall, and therefore we must have g(r) — > 
as r — > r s . Using equation (47) to substitute for v 2 , we can linearize the factor in parentheses in 
equation (46) to obtain 

■9vl £ 2 Q 7 ( r 2 v 2 s \^ /I + 7 . \ (m 
Solving this equation for g(r) yields the asymptotic result 



2e r 2 -£ 2 -( 1 + l)f{r 

y v ) - 

where 



r 2 v£(r) - (7 - l)/(r) 



_ 2 7 r 2 
/(r) = 



K 2 



2r 3 (r — r s ) 



7-1 



7+1 

(50) 
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Note that the function f(r) diverges as r — ► r s , but it does so much more slowly than Vg(r), and 
consequently g(r) — > at the horizon as required. Equations (47), (49), and (50) provide a useful 
asymptotic representation for the inflow velocity v(r) that describes the first-order correction to 
purely free-fall behavior close to the horizon. 

By combining equations (8), (45), and (47), we can show that the asymptotic solution for a{r) 
is given by 



a{r) = 

where 



^r) [l + «W]'(r-r a r, (51) 

"wn-)- <52) 

The dominant behavior as r — > r s is a oc (r — r s ) CT . The exponent a is negative, and therefore a 
diverges at the horizon, albeit much more slowly than v, which approaches v g oc (r — r s ) -1 / 2 . 

We can also easily determine the leading behavior of the disk half-thickness H(r) close to the 
horizon by combining the hydrostatic relation H = 6o<V^K with equations (8) and (51), which 
yields 



where c is the speed of light and 

As r — > r s , the dominant behavior is i7(r) oc (r — r s ) 5 . Since the flow becomes adiabatic close to 
the horizon, we can use equations (7) and (51) to show that the dominant asymptotic forms for 
the pressure P and the mass density p are given by 

P(r) oc (r — r s ) A , p(r) oc (r — Tg)' 7 , (55) 

where 

A = - -1- , r? = - — — . (56) 
7+1 7+1 

Table 1 includes values for the exponents a, 5, A, and rj obtained for several different values of the 
adiabatic index 7. Note that 5 ~ 1, indicating that iT is roughly proportional to r — r g . Hence the 
disk has zero thickness (i.e., a cusp) at r = r s , reflecting the fact that the gas pressure is unable to 
support the disk against the strong gradient of the gravitational potential as the matter approaches 
the horizon. This particular behavior is a manifestation of the assumption of vertical hydrostatic 
equilibrium. In § 6.2 we explore the consequences of replacing this assumption with radial free-fall 
close to the horizon. 

Equations (47) and (51) provide extremely accurate asymptotic solutions for v(r) and a(r), 
respectively. The accuracy of these solutions will be demonstrated in § 5 by comparing them with 
exact solutions obtained by integrating numerically the differential conservation equations. Taken 
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together, equations (38), (43), (47), and (51) completely determine the nature of the flow close to 
the horizon in terms of the three free parameters £q, Kq, and eo, which describe the specific angular 
momentum, the specific entropy, and the specific energy of the particles entering the black hole, 
respectively. As discussed in § 5, the requirement of smooth passage through a critical point imposes 
an additional constraint that effectively reduces the number of free parameters from three to two. 
The asymptotic results we have derived in this section can be used to define boundary conditions 
applicable close to the event horizon that serve as the basis for numerical simulations of the global 
structure of advection-dominated disks. In § 5 we perform global disk structure calculations and 
compare the numerical solutions obtained with our asymptotic expressions in the vicinity of the 
horizon. 



The various differential and algebraic conservation equations may be solved numerically to 
determine the profiles of v, a, and t = r 2 Q. Several methods are available to solve the coupled 
system of equations, such as explicit integration using a Runga-Kutta solver, or the utilization of a 
global, iterative relaxation method (e.g., Press et al. 1986). Each of these techniques requires the 
imposition of boundary conditions at the edges of the computational domain. The new asymptotic 
relations derived in § 4 can be applied at a radius just outside the horizon to provide the inner 
boundary conditions at the black hole event horizon needed to compute global solutions for the disk 
structure. Before attempting to solve the computational problem to determine the disk properties, 
it is worthwhile to review the critical nature of the conservation equations. 



Black hole accretion flows are in general transonic, and therefore the computational domain 
can be broken into two regions, one above the sonic radius and one below it. Successful global 
solutions must pass smoothly through the sonic point, which is a critical point for the flow. In 
order to explore the critical nature of the flow, it is convenient to derive a dynamical equation 
based on the mass, momentum, and energy conservation equations. As a preliminary step, we can 
use equations (3), (6), and (7) to express the density p in terms of v, r, and P as 



Using this relation to eliminate p in the entropy equation (14), we can derive an equation for the 
pressure derivative dP/dr. The result obtained is 



5. 



GLOBAL FLOW SOLUTIONS 



5.1. Dynamical Equation and Critical Conditions 




(57) 




(58) 
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This can be used to eliminate the pressure derivative in the radial momentum equation (9) to yield 
the dynamical equation 

(v 2 2 1 \d\nv_ l 2 -e K 2 7 /3 d\n£ K \ ( 1 _ - 1 \ v £ K {£ - £ ) 2 
\a 2 7+iy dr a 2 r 3 7 +l\r dr J V7+I/ «« 4? " 4 

where 



(GM) 1 / 2 r 3 / 2 

£ K (r) = r 2 n K (r) = (60) 

r — r s 

denotes the specific angular momentum of particles in circular, Keplerian orbits at radius r, and 
we have also used equations (8) and (29). Equation (59) agrees with equation (2.16) of Narayan 
et al. (1997). Global flow solutions can be obtained by integrating simultaneously the two coupled 
differential equations (30) and (59), which govern the functions £{r) and v(r), respectively. This is 
similar to the procedure followed by Narayan et al. (1997), except that they included an additional 
differential equation for a(r), based on the entropy equation (14). However, this extra differential 
equation is not necessary because the energy flow rate E is conserved when radiative losses are 
neglected, as assumed in the ADAF scenario (Molteni, Gerardi, & Valenza 2001). This fact allows 
us to solve for a as an algebraic function of v, £, and r using equation (44), which yields 



7 -1 / _ e_ _ £o£ , GM 

^r\ €o ~Y + 2^ 2 ~^ 2 1 



« 2 = 1 Uo-Tr + ^-^T + • (61) 



Critical points occur where the numerator and denominator in equation (59) for dlnv/dr 
vanish simultaneously. This yields the critical conditions 

v 2 2^ 

- - — — r = , r = r Cl (62) 
a z 7 + I 



a 2 r 3 



• 3 7+ 1 \r dr J V7+1/ aa r 



where r c is the critical radius. Global solutions must pass through a critical point, and therefore 
equations (61), (62), and (63) can be used to interrelate the six quantities (r c , v c , a c , £ c , eo, £0), where 
v c , a c , and £ c denote quantities measured at the critical radius r = r c . We can integrate the system 
of equations (30) and (59) away from the critical point, either towards large radii or towards the 
horizon. However, due to the nature of the critical point, we cannot begin the integration precisely 
at r = r c . We must therefore employ l'Hopital's rule to evaluate dv/dr at the critical point, and 
then perform a linear extrapolation to offset the starting conditions slightly in radius (Molteni et al. 
2001; Chen & Taam 1993; Chen et al. 1997). This procedure involves the solution of a quadratic 
equation for the critical value of dv/dr. In our application, the negative, real root gives the value 
for the derivative at the critical point. 
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5.2. Disk Structure Calculations 

In order to illustrate the utility of the asymptotic relations developed in § 4, we shall per- 
form several calculations of the disk structure based upon explicit integration of the differential 
equations. Since there are two coupled differential equations in the system, there are two linearly 
independent local solutions around the singular point at the horizon. Only one of the local solutions 
is physically acceptable, and this is the solution that we have obtained asymptotic representations 
for in § 4. Explicit integration of the system of equations from the critical point towards the horizon 
is unstable, because in general it is impossible to avoid exciting the second, linearly independent 
solution, which displays an unphysical behavior as the gas approaches the horizon. While it is an 
open question whether shocks supported by a "centrifugal barrier" exist in black hole accretion 
disks (e.g., Chakrabarti 1997), our goal here will be to develop global, shock- free solutions in order 
to demonstrate the utility of the boundary conditions derived in § 4 in the simplest possible manner. 

With the availability of the asymptotic expressions for v(r) and t(r) derived in § 4, we can 
employ an explicit, stable integration in the outward direction, starting at a point just outside 
r = r s . In this approach, equations (38) and (47) are used to set the starting values for £ and v, 
respectively, as functions of the three constants (eo, £o, Kq). Numerical integration of equations (30), 
(59), and (61) in the outward direction yields solutions for £(r), v(r), and a{r). For given values of 
£o and eo, the parameter Kq can be determined by requiring that the flow pass smoothly through 
a critical point, at some radius r = r c . In order to determine £q and eo, we must therefore supply 
two additional boundary conditions. These extra conditions are usually imposed by requiring that 
the disk become Keplerian and geometrically thin at some arbitrary outer radius (Narayan et al. 
1997). However, it is not completely clear whether ADAF solutions can merge smoothly with cool 
thin disks (Yuan 1999; Yuan et al. 2000). This particular issue is not central to our considerations 
in this paper, since our focus here is on discussing the inner boundary conditions appropriate for 
advection-dominated black hole accretion. Therefore, in order to develop numerical examples that 
illustrate the role of the inner boundary condition without undue complexity, we shall simply set 
eo = and require that a — > as r — > oo, in keeping with the self-similar, advection-dominated 
models (Narayan & Yi 1994, 1995; Blandford & Begelman 1999; Becker, Subramanian, & Kazanas 
2001). The self-consistency of the numerical solutions obtained for £(r), v(r), and a(r) using the 
outward integration is checked by confirming that at the critical radius, £(r c ) = £ c , v(r c ) = v c , and 
a(r c ) = a c , where the quantities (v c , a c , £ c , r c ) satisfy the critical conditions given by equations (61), 
(62), and (63). 

The integrations begin at a starting radius, r*, located close to the event horizon, and proceed 
in the outward direction, towards the critical point. In our numerical examples, we shall work in 
terms of natural gravitational units (GM = c = 1, r s = 2), and the starting radius will be given 
by r* = 2.001, which is just outside the Schwarzschild radius. The corresponding starting values 
for the specific angular momentum £* = £(r*) and the inflow velocity = u(r*) are computed by 
applying the asymptotic formulas given by equation (38) and (47), respectively, at the radius r = r*. 
The starting value for the isothermal sound speed, a* = a(r*), is determined using equation (61). 
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The exact solutions for v, a, and I as functions of radius are obtained by integrating numerically 
the coupled equations (30), (59), and (61). The values of the various model parameters are listed 
in Table 2 for the three models that we consider in detail below. The specific parameter values we 
have selected correspond to those adopted by Narayan et al. (1997) in several of their calculations. 

In Figures 1, 2, 3, 4, and 5, we present solutions obtained by setting £q = 2.6, eo = 0, 
a = 0.1, and 7 = 1.5, which we refer to as Model 1. Note that this value for 7 represents 
approximate equipartition between the gas pressure and the magnetic pressure. In this scenario, 
the a-prescription we have employed for the viscosity (eq. [11]) is valid even within the supersonic 
region of the flow. The exact numerical solution for the inflow velocity v (r) is compared with 
the asymptotic solution given by equation (47) in Figure 1. Also included for comparison are 
the free-fall velocity distribution, vs(r) = [2GM/(r — rg)] 1 / 2 , and the exact numerical solution 
for the isothermal sound speed a(r). The agreement between the exact solution for v(r) and the 
asymptotic expression is excellent for 2 < r < 3, and remains reasonably close all the way out 
to the critical point, which is located at r c = 6.132 for this model. Note that v remains well 
below Vfi until the material gets quite close to the horizon. The exact numerical solution for the 
sound speed a(r) is compared with the asymptotic solution (eq. [51]) in Figure 2. The agreement 
between these two results is surprisingly close all the way out to the critical radius. In Figure 3 
we plot the exact numerical solution for the specific angular momentum, I = r 2 fi, along with the 
asymptotic formula given by equation (38). The two expressions for £{r) merge smoothly asr —> r s , 
in validation of our asymptotic analysis. In Figure 4 we display the global solutions obtained by 
joining the numerical results for v(r) and a(r) from Figure 1 with curves generated by integrating 
in the outward direction starting from the critical radius r = r c . Note that the global solutions pass 
smoothly through the critical point as required. In Figure 5 we display the results obtained for v(r) 
and a(r) by integrating in the inward direction starting from the critical point. These are plotted 
along with the numerical solutions obtained using the outwardly-directed integration starting close 
to the horizon (Fig. 1). Near the critical point, the two sets of solutions agree closely. However, 
the inwardly-directed integration is unstable, and the denominator of the dynamical equation (59) 
vanishes at r ~ 5, where v 2 = 27 a 2 /(l +7). The existence of this instability provides one of the 
main motivations for developing the inner boundary conditions in § 4, since the availability of these 
boundary conditions facilitates the integration in the outward direction, which is stable. 

In Figures 6, 7, 8, and 9, we present results for Model 2, with £q = 1.76, eo = 0, a = 0.3, 
and 7 = 1.5. The exact and asymptotic solutions for the inflow velocity v(r) are compared in 
Figure 6. Note that two results are indistinguishable in the entire region between the horizon 
and the sonic (critical) point, which is located at r c = 10.63 for the Model 2 parameters. The 
exact solution for the sound speed a(r) obtained via numerical integration is compared with the 
asymptotic expression (eq. [51]) in Figure 7. In this case, the two results agree closely for 2 < r < 3, 
although the agreement deteriorates for r ^ 5. The exact numerical solution for the specific angular 
momentum t = r 2 Q, is plotted in Figure 8 along with the asymptotic result given by equation (38). 
The two expressions again display a smooth merger as the gas approaches the horizon. In Figure 9 
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we plot complete global solutions obtained by combining the Figure 6 results with solutions for 
v(r) and a{r) obtained by integrating away from the critical point in the outward direction. The 
global solutions pass smoothly through the critical point, in satisfaction of the critical conditions. 

It is interesting to compare our results with those obtained by Narayan et al. (1997). An 
examination of their Figures 1, 2, and 3 reveals close agreement with our results. In fact, their 
plots of the variation of £ near the horizon all show that dijdr — > as r — > r s , which is consistent 
with our prediction based on equation (38). This is intriguing considering the fact that they did 
not formally adopt boundary conditions identical to ours in their calculations. In fact, they state 
that d£l/dr = at the horizon in the discussion following their equation (2.18), which is obviously 
incorrect. Hence their treatment of the inner boundary conditions is unclear. Our numerical results 
also agree with those of Chen et al. (1997), although these authors do not extend their calculations 
to the horizon, and instead truncate the disk at an arbitrary radius r- m = 3. For the starting radius 
in our outwardly-directed integrations, we have set r* = 2.001, and the corresponding local free-fall 
velocity there is vg(r*) = 44.7214. The values we obtain for the starting velocity displayed 
in Table 2 are just slightly smaller than the local free-fall velocity, as a result of centripetal and 
pressure effects. However, this small difference is crucial for determining the velocity distribution 
via the subsequent integration away from the horizon. 

Along with the results associated with Models 1 and 2 results, Table 2 also includes a summary 
of the results obtained for Model 3, with £q = 3.21, eo = 0, a = 0.03, and 7 = 1.5. It is interesting 
to compare the values obtained for the entropy and angular momentum parameters in the three 
models. In particular, Table 2 includes results for the entropy function K(r) (eq. [15]) evaluated 
at the critical radius, K c = K(r c ), as well as at the horizon, Kq = K(r s ). Note that Kq always 
exceeds K c , reflecting the fact that viscous dissipation between the sonic point and the horizon 
increases the entropy of the gas. Our results indicate that Kq exceeds K c by 91% for a = 0.3, by 
37% for a = 0.1, and by 14% for a = 0.03. This is consistent with the expectation that larger 
values of a should be associated with higher dissipation rates. Also note that the accreted specific 
angular momentum, £q, is always lower than the specific angular momentum at the sonic point, 
£ c , and that the ratio £q/£ c decreases with increasing a, as expected, due to the increase in the 
magnitude of the viscous torque. 

6. DISCUSSION 

The development of global solutions for the structure of advection-dominated accretion disks 
requires careful consideration of the boundary conditions that apply close to the black hole event 
horizon, because the horizon always represents the fundamental inner boundary for any such cal- 
culation. In this paper we have explored the consequences of the pseudo-Newtonian potential for 
the asymptotic structure (in the vicinity of the horizon) of advection-dominated accretion disks 
incorporating the Shakura-Sunyaev viscosity prescription. The presence of the pseudo-Newtonian 
potential introduces a regular singular point into the differential conservation equations. In recog- 
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nition of this fact, we have employed asymptotic analysis to determine the leading behaviors of the 
physical quantities v, a, £, and K close to the horizon. The main asymptotic results derived in 
§ 4 are given by equation (38) for £(r), equation (43) for K(r), equation (47) for v(r), and equa- 
tion (51) for a{r). These expressions clearly illustrate how the various physical quantities depend 
on the three parameters eo, £q, and K , denoting the accreted specific energy, the accreted specific 
angular momentum, and the accreted specific entropy, respectively. 

The analytical, asymptotic relations we have obtained can be used to provide the inner bound- 
ary conditions needed for simulation of the disk structure. These boundary conditions are essential 
whether the computation is performed using an explicit integration or an iterative relaxation tech- 
nique. In § 5.2, the new boundary conditions were used to perform an explicit numerical integration 
of the coupled conservation equations in the outward direction, starting at a point just outside the 
event horizon and ending at the critical point. The exact numerical solutions obtained for v, a, 
and i were compared with the corresponding analytical expressions derived in § 4. The agreement 
between the two sets of solutions is remarkable, providing positive confirmation of the validity of 
our asymptotic analysis and the resulting boundary conditions. We have demonstrated in Ap- 
pendix A that the velocity distribution associated with our pseudo-Newtonian model agrees with 
general relativity in the vicinity of the horizon. In the remainder of this section, we will examine 
a few additional questions related to the self-consistency of the model. 



6.1. Torque Boundary Condition 

In our model, the torque Q vanishes at the horizon, rather than at the radius of marginal 
stability, r ms = 6 GM/c 2 . This is consistent with simulations performed by Hawley & Krolik (2001), 
Agol & Krolik (2000), and Gammie (1999), who all find that the stress has a finite value at r = r ms . 
Furthermore, our results demonstrate that any advection-dominated disk with hydrostatic vertical 
structure must have di/dr = at the horizon. This in turn implies that the radial derivative of 
the torque must vanish there, i.e., dQ j dr = at r = r s (eq. [41]). This new condition supplements 
the well-known requirement that Q = on the horizon. We show below that this behavior can be 
understood as a simple consequence of the adiabatic nature of the flow close to the horizon. First 
we combine the fundamental expression for the torque, equation (10), with the mass conservation 
equation (3) and the Shakura-Sunyaev viscosity prescription (eq. [11]) to obtain 

G = - aMa2r ^ . (64) 

v n K dr y 1 

We have found that close to the horizon, the viscous heating vanishes and therefore the pressure 
obeys the adiabatic relation P oc p 7 , implying that a? oc p 7_1 . Substituting for p using the mass 
conservation relation M = AirrHpv, we find that a 2 oc (rHv) 1 ^ asr^r s . In a hydrostatic disk, 
H oc a/f^K, and consequently a oc (rv/ft K ) (1 ~ 7)/{1+7) (cf. eq. [45]). Since d£l/dr approaches a 
finite value as r — > r s , and v approaches free-fall, we immediately conclude that Q oc (r — r s )^, with 
(3 = (7 + 5)/ (2 + 27). This agrees with our earlier derivation (see Table 1), and confirms that (3 > 1, 
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with the resulting implication that dQ jdr must vanish at the horizon. Hence we have demonstrated 
that the vanishing of the derivative of the torque at the horizon is a predictable consequence of the 
scalings of a and v as r — > r s in standard ADAF disks. 

6.2. Effects of Central Free-Fall 

In this paper we have focused on the behavior of advection-dominated disks that are in ver- 
tical hydrostatic equilibrium at all radii, because we are interested in developing inner boundary 
conditions that are consistent with the standard ADAF scenario (Narayan et al. 1997; Chen et al. 
1997). The boundary conditions and asymptotic behaviors we have obtained agree well with global 
ADAF simulations, and therefore they provide a useful foundation for numerical calculations of the 
disk structure. While our approach is a reasonable strategy from a computational point of view, 
we should acknowledge that in reality, the gas will probably stop responding to vertical pressure 
forces in the supersonic region close to the horizon. It is perhaps more likely that the gas crosses 
the horizon in nearly radial free-fall, with H oc r. From a physical point of view, it is worthwhile 
to consider how this alternative central inflow condition would affect some of the basic conclusions 
we have reached in this paper. 

We have shown quite generally in § 4.3 that the dissipation rate vanishes as the gas approaches 
the horizon, and therefore P oc /? 7 and a? oc p 7_1 . This result does not depend on the assumption 
of vertical hydrostatic equilibrium. Proceeding as in § 6.1, we substitute for p using the mass 
conservation equation (3) and find that, in general, a? oc (rHv) 1 ^' 1 as r — > r g . In the quasi- 
spherical free-fall region close to the horizon, the disk half-thickness is given approximately by 

H = d r , (65) 

where do is a constant. This implies that asr^r s , the variation of the sound speed satisfies 

a 2 oc {r 2 v) 1 ^ . (66) 

In the free-fall case, it is convenient to define the entropy function using the alternative form 

K( r ) = r 2 ua 2/( 7 -l) _ (g7) 

Although differs from the definition of K(r) used in the hydrostatic scenario (eq. [15]), by combining 
equations (3), (7), (65), and (67), we can confirm that K 1 ^ 1 oc Up^' y , and therefore K is a linear 
function of the entropy per particle S if the disk is in free-fall (see eqs. [16] and [17]). Next we 
express the asymptotic behavior of £(r) close to the horizon using 

i{r) = £ + B(r-r s f , (68) 

where £q is the accreted specific angular momentum, and the constants B and f3 are analogous to 
the hydrostatic constants B and (3 appearing in equation (31). Following the same steps as in § 4.3, 
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we now obtain 

Using asymptotic analysis and requiring that exponents and coefficients balance on the two sides 
of this expression, it is straightforward to show that the solutions for B and (3 are given by 



where c is the speed of light and 



K = lim K{r) (71) 



denotes the value of the entropy function at the horizon. Note that f3 > 1 for all 4/3 < 7 < 5/3, 
and consequently dQ/dr = di/dr = at the horizon. This is the same conclusion reached in § 4.3 
under the assumption of vertical hydrostatic equilibrium, although the exponent (3 differs slightly 
from the hydrostatic exponent (3 = (7 + 5)/(2 + 27). Hence the vanishing of the derivative of the 
torque at the horizon in ADAF disks is a very general result. 



6.3. Causality of Viscous Transport at the Horizon 

Numerous authors have pointed out that the diffusive nature of the angular momentum trans- 
port associated with the Shakura-Sunyaev viscosity prescription v = aa 2 /^K can lead to causality 
violations in accretion disks (e.g., Kato 1994; Narayan 1992). This issue can be most easily under- 
stood by considering the evolution of an initially localized component of the angular momentum 
distribution, represented by a 5- function at some arbitrary radius r = vq and arbitrary time t = to- 
As time proceeds, the distribution will spread in radius in an approximately Gaussian manner, im- 
plying propagation to infinite distance in a finite time, which violates causality. This phenomenon 
has a negligible effect on the structure of the disk in the outer, subsonic region, because the mean 
transport velocity for the angular momentum is typically very small despite the fact that an in- 
finitesimal portion of the signal propagates with infinite speed. However, the question of causality 
needs to be examined carefully in the inner, supersonic region, where the radial inflow velocity v r 
approaches the speed of light, and all signals should therefore be advected into the black hole. In 
order to address this issue in the context of the ADAF scenario considered here, it is useful to ex- 
amine the diffusion equation governing the angular momentum distribution in the disk. Following 
Blandford Sz Begelman (1999), we write the time-dependent equation as 

5^-^-0, (72) 

where 

u = -iirrp = — (73) 
v 
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represents the mass per unit radius in the disk. By combining equation (72) with equation (10) for 
the torque Q, we can obtain the alternative form 



dL_d_ 
dt dr 



dL vL dn 2vL 
vL + u- -f 

or fi or r 



(74) 



where 

L = fir 2 n (75) 
denotes the angular momentum per unit radius. 

To obtain further insight, we can recast equation (74) in the form of the Fokker-Planck equation 

91 - a WL) + *(l"£), (76) 



dt dr \ dt J dr 2 \ 2 dt 
where the Fokker-Planck coefficients 

d(r) 2v dv v d\x 1 do 2 

nr = T + d-r + ]idr-- v > 2""' (77) 

describe, respectively, the rates of "drifting" and "broadening" experienced by the initially localized 
angular momentum distribution due to diffusion (Reif 1965). We are interested in evaluating the 
Fokker-Planck coefficients in the context of the steady-state ADAF disks considered here, which 
have jU oc v^ 1 . In the inner region, close to the horizon, we have found that v oc vq = [2GM/(r — 
r-g)] 1 / 2 and a oc (r — r s )( 1-7 )/( 27+2 ). Consequently, the Shakura-Sunyaev viscosity, v = aa 2 /^}^, 
displays the asymptotic behavior v oc (r — r s ) 2 ^ 1+ ^ in the vicinity of the horizon. Since v vanishes 
as r — > r s , it follows that do 2 /dt — > 0. One can also show that the mean transport velocity, 
d(r)/dt, approaches —v as r — > r s . The vanishing of the "broadening" rate at the horizon implies 
that the angular momentum is simply advected into the black hole, and there is no nonphysical 
transport to infinite distance in finite time. Taken together, our asymptotic results for the Fokker- 
Planck coefficients da 2 /dt and d(r)/dt demonstrate that the transport of angular momentum at 
the horizon is causal, in agreement with general relativity. We have therefore confirmed that, 
at least in the context of the ADAF model considered here, there are no causality violations at 
the horizon associated with the Shakura-Sunyaev prescription for the viscosity. Interestingly, this 
result remains valid even if the hydrostatic relation is replaced with the central free-fall condition 
discussed in S 6.2. 



6.4. Energy and Mass Loss 

In this paper, we have adopted the perfect ADAF approximation by assuming that all of the 
transfer rates E, J, and M are constant. We have therefore neglected the possibility of radiative 
losses, as well as outflows of matter and energy. This contradicts recent observations suggesting 
that many low-luminosity X-ray AGNs possess relativistic jets that originate very close to the 
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event horizon of the central black hole (Fender 2001; Nagar et al. 2002). Examples include RS 
1915+105 (Belloni et al. 1997; Dhawan, Mirabel, k Rodriguez 2000), XTE J1118+480 (Fender et 
al. 2001), XTE J1550-564 (Corbel et al. 2001), GS 1354-64 (Brocksopp et al. 2001), and perhaps 
Cyg X-l (Stirling et al. 2001). The fact that disks and jets often appear together suggests that 
the presence of the outflows may be a necessary ingredient for the accretion to proceed. This 
possibility has motivated several investigations into the relationship between outflows and ADAFs. 
From a theoretical viewpoint, the positivity of the Bernoulli parameter in advection-dominated flows 
suggests that the gas in these systems is gravitationally unbound (Narayan, Kato, & Honma 1997; 
Narayan & Yi 1994, 1995). Based on this observation, Blandford & Begelman (1999) investigated 
the effects of mass, energy, and angular momentum loss on the structure of the disk in the context of 
a self-similar model incorporating Newtonian gravity. This approach was extended by Becker et al. 
(2001) to describe self-similar disk/outflow systems governed by the pseudo-Newtonian potential. 

Observations of strong radio emission from X-ray underluminous AGNs suggest that relativis- 
tic particles abound in the hot plasmas. The typical energy of these particles is much higher than 
the average thermal energy of the gas, implying the presence of an efficient acceleration mecha- 
nism. In this connection, it is interesting to note that the low density in advection-dominated 
disks makes them plausible sites for the acceleration of relativistic particles via interactions with 
magnetohydrodynamical waves because the plasma is so tenuous that ion-ion collisions are unable 
to thermalize the energy of the accelerated particles (Becker et al. 2001; Subramanian, Becker, & 
Kazanas 1999). Hence in the X-ray underluminous AGNs, particle acceleration and the resulting 
outflows of unbound particles from the disk may represent the dominant cooling mechanism, re- 
moving excess energy and thereby allowing accretion to proceed. Conversely, in the X-ray bright 
systems, the efficiency of particle acceleration in the disk is lower due to the higher density, which 
tends to thermalize the energy of the accelerated particles. In these systems, it is the X-ray emis- 
sion that removes most of the binding energy. This interpretation helps to explain the observed 
anticorrelation between the outflow strength and the X-ray luminosity, as well as the positive cor- 
relation between the X-ray hardness ratio and the radio emission strength (Celotti & Blandford 
2001; Corbel et al. 2000). 

While outflows have not been incorporated into our analysis, we expect that their inclusion 
would have little if any effect on the boundary conditions at the event horizon that we have derived. 
This is because the power source for the outflows would presumably be the viscous dissipation, which 
clearly vanishes rather quickly below r = r ms . Hence the physics in the asymptotic region close 
to the horizon should be insensitive to the production of the outflows. In future work, we plan 
to utilize our asymptotic relations to facilitate the development of detailed disk models including 
outflows of matter and energy that are self-consistently coupled with the disk. 
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6.5. Conclusion 

In this paper we have obtained a number of useful analytical expressions that completely 
describe the structure (close to the event horizon) of advection-dominated, pseudo-Newtonian ac- 
cretion flows based on the Shakura-Sunyaev viscosity prescription. The dynamical results depend 
on three quantities: the accreted specific energy eo, the accreted specific angular momentum £q, 
and the accreted specific entropy klnK®. In our approach, eo and Iq are treated as free parameters, 
and the value of Kq is determined by requiring that the global solution pass smoothly through a 
critical point. The asymptotic expressions derived in § 4 provide a set of inner boundary conditions 
that can serve as the basis for subsequent numerical integration of the conservation equations. 
We emphasize that any physically consistent one-dimensional advection-dominated accretion disk 
model based on the Shakura-Sunyaev viscosity prescription must satisfy these boundary conditions. 
The analytical expressions agree extremely well with the exact numerical solutions out to a few 
gravitational radii from the horizon. Hence our results provide a valid description of the essential 
physics of the accretion process in the vicinity of the horizon. 

We have assumed that the effects of general relativity can be approximated using the pseudo- 
Newtonian gravitational potential. While this is a widely used approximation that preserves many 
of the important dynamical characteristics of flows in the Schwarzschild metric, one may well ask 
how the specific results we have obtained here translate into full general relativity. For example, 
will the radial derivative of the torque really vanish at the horizon in the Schwarzschild metric? 
Obviously this question cannot be answered definitively without employing a fully relativistic cal- 
culation. However, we demonstrate in Appendix A that the motions of particles near the horizon 
predicted by our pseudo-Newtonian model are in complete agreement with the actual motions of 
freely-falling particles in the Schwarzschild metric. Furthermore, we have established in § 6.3 that 
the diffusive transport of angular momentum in our model is causal in nature at the horizon. This 
confirms that pseudo-Newtonian ADAF models incorporating the Shakura-Sunyaev viscosity pre- 
scription can be used to describe the structure of the disk all the way to the event horizon. We 
argue based on these results that the general characteristics of the asymptotic solutions we have 
obtained are likely to be preserved in full general relativity. 

The numerical examples presented in § 5 do not include shocks, which may occur in accretion 
flows as a result of the interaction between the gas and a "centrifugal barrier" (Chakrabarti 1997). 
Shocks may also play a role in powering the outflows associated observationally with hot disks 
(Yuan et al. 2002). Although shocks were not explicitly considered in our analysis, we argue 
that our asymptotic results should apply equally well whether or not shocks are present, because 
our results are based on the fundamental physical processes operative near the horizon, and those 
processes are insensitive to the history of the gas. A related question concerns the relevance of our 
results in the context of convection-dominated accretion flows ("CDAFs"). These are stationary, 
convective-envelope solutions that technically have zero accretion rates, although in fact a small 
amount of matter is expected to flow into the black hole (Narayan, Igumenshchev, & Abramowicz 
2000). We emphasize that our basic results should apply in this situation as well because the gas 
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that enters the black hole must satisfy the same asymptotic conservation equations, independent 
of whether it has passed through a CDAF, a shocked disk, or a conventional ADAF. 

The authors are grateful to the anonymous referee for helping us to improve the discussion. 
PAB would also like to acknowledge several stimulating conversations with Menas Kafatos, Ken 
Wolfram, and Demos Kazanas, as well as generous support from the Naval Research Laboratory 
during a portion of the research. 
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A. APPENDIX 

PSEUDO-NEWTONIAN PARTICLE DYNAMICS 

In our approach to modeling the disk structure, we have incorporated the effects of general 
relativity in an approximate manner by utilizing the pseudo-Newtonian potential, 



$(r) 



-GM 



r — r a 



(Al) 



which gives correct results for the radius of marginal stability (r ms = 6GM/c 2 ), the marginally 
bound radius (r mD = 4GM/c 2 ), and the horizon radius (r s = 2GM/c 2 ) around a nonrotating 
black hole (Paczyhski & Wiita 1980; Abramowicz, Calvani, & Nobili 1980). In order to understand 
how the motions of particles in the pseudo-Newtonian potential are related to the exact solutions 
given by general relativity, we shall briefly review the dynamics of particles freely falling in the 
Schwarzschild metric. We begin by writing down exact expressions for the radial velocity v r and 
the azimuthal velocity describing the motion of a particle as measured by a local, static observer 
in the Schwarzschild geometry. Using equations (12.4.17) and (12.4.18) from Shapiro & Teukolsky 
(1983), we obtain 



and 



/ r s \V2 £c 



(A3) 



where E and i denote the particle's specific energy and specific angular momentum at infinity, 
respectively. The locally measured specific energy, E\ oca y, is related to E by 



-EWai = e{\ — -^j 
In terms of the locally measured Lorentz factor, 



'1/2 



Tlocal — 



-Elocal 



(A4) 



(A5) 



the radial component of the four-velocity, v r , measured by a static observer at radius r can be 
written as 



_ r r E ( r s \-i/2 

V r = Tlocal V = — 1 - 

c V r 



Z7 \ -2 

E\ f l _r ± 

r 



1 + 



c 2 r 2 



nl/2 



and the locally measured azimuthal component of the four- velocity, v v , is given by 

e 



V v = ri ocal V 



¥ — 



(A6) 



(A7) 
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We shall now focus on the case of a particle falling from rest at infinity, with E = c 2 . Note 
that I can still have an arbitrary value in this case, since the azimuthal velocity vanishes at infinity 
for all values of £. Our expression for v r now reduces to 

2 2 GM £ 2 

V 2 r = 2 , A8 

r — r s r z 

or, equivalently, 

\v 2 r + \vl + ${r) = 0, (A9) 

where $(r) is the pseudo-Newtonian potential given by equation (Al). Equation (A9) resembles a 
classical Newtonian energy equation, except that v r and v v are four-velocities rather than conven- 
tional velocities. This is, in fact, one of the basic motivations for introducing the pseudo-Newtonian 
potential. Despite its classical appearance, the result obtained for v r by solving equation (A9) with 
v v = £/r is exactly equal to the radial four- velocity of a particle freely falling from rest at infinity 
in the Schwarzschild metric. Note that as the particle approaches the event horizon of the black 
hole, the radial four- velocity diverges, i.e., 

2GM\ 1/2 



(2GM\ 
\r-rj 



v r (r) -»■ v s (r) = , r -> r s . (A10) 



On the other hand, the azimuthal four- velocity, v v , remains bounded, and we find that for a given 
value of £, 

£ 

Vy{r) -> — , r -> r s . (All) 

r s 

As the particle approaches the horizon, the limiting values for the physical velocities v r and are 
given by 

lim v f = c , lim = , (A12) 

r— »r g r ^ r s 

which follow from equations (A2) and (A3). Hence a stationary observer close to the horizon sees 
the particle falling radially inward at the speed of light, as expected. 

As a consequence of utilizing the pseudo-Newtonian potential to describe ADAF disks, we have 
found in § 4 that the asymptotic behavior of the "radial velocity" v is given by 

'2GM\ 1/2 

— r > r^r s , (A13) 
' s J 

and the asymptotic behavior of the "azimuthal velocity" w = r f2 is likewise given by 

w -> ^ , r^r s , (A14) 

where £q = J/M is the specific angular momentum of the material crossing the event horizon. 
Comparing these expressions with equations (A10) and (All), we find that close to the horizon, the 
quantities v and w, respectively, are exactly equal to the radial component (v r ) and the azimuthal 
component (v^) of the four- velocity for a particle freely falling in the Schwarzschild metric. 
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Fig. 1. — Exact numerical solution for the inflow velocity v (solid line) plotted as a function of 
radius r for Model 1, with a = 0.1 and £q = 2.6. Included for comparison is the asymptotic 
solution for v(r) given by eq. (47) (dashed line). The two results agree closely for 2 < r < 3, 
and remain similar out to the critical point at r c = 6.132. Also shown are the free-fall velocity 
vg = [2GM/(r — Tg)] 1 / 2 (dotted line), and the numerical result for the sound speed a multiplied by 

[27/(7 + l)] 1 ^ 2 (dot-dashed line). The sound speed curve crosses the numerical solution for v at the 
critical point. 
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Fig. 2. — Exact numerical solution for the isothermal sound speed a (solid line) plotted as a function 
of radius r for Model 1. Parameters for this model are listed in Table 2. Included for comparison 
is the analytical, asymptotic solution for a(r) (eq. [51]; dashed line). The two results agree well 
between the starting radius r* = 2.001 and the critical point r c = 6.132. 
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Fig. 3. — Exact numerical solution for the specific angular momentum i = r 2 O plotted as a function 
of radius r (solid line) for Model 1. For comparison we include the asymptotic solution for £(r) 
given by eq. (38) (dashed line). The functions merge smoothly as r — > r s , and approach the accreted 
specific angular momentum, £q = 2.6 (dotted line). 
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Fig. 4. — Global numerical solutions for the inflow velocity v(r) (solid line) and the sound speed 
a multiplied by [27/(7 + l)] 1 ^ 2 (dashed line) plotted as functions of radius r for Model 1. The v 
and a curves pass smoothly through the critical point at r c = 6.132, where they cross. The flow is 
supersonic for r < r c . The vertical dotted line indicates the location of the critical point. 
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Fig. 5. — Numerical solutions for the inflow velocity v(r) (solid line) and the sound speed a 
multiplied by [27/(7 + l)] 1 / 2 (dashed line) plotted as functions of radius r for Model 1. Separate 
results are indicated for the inwardly and outwardly directed integrations. Note that the two sets 
of results agree near the critical radius at r c = 6.132. However, the inwardly-directed integration 
fails at r ~ 5, where v 2 = 2 7 a 2 /(I + 7). See the discussion in the text. 
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Fig. 6. — Same as Fig. 1, except results correspond to Model 2, with a = 0.3, lo = 1.76, and 
r c = 10.63. Note that in this case, the asymptotic solution for v(r) is indistinguishable from the 
exact numerical solution in the entire region between the horizon and the critical point. Additional 
parameters for this model are listed in Table 2. 
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Fig. 7. — Comparison of solutions for the isothermal sound speed a(r) obtained in Model 2. In this 
case, the exact numerical solution (solid line) agrees with the analytical solution (dashed line) close 
to the horizon, but diverges for r > 5. The location of the critical radius at r c = 10.63 is indicated. 
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Table 1. Exponents of Solutions as Functions of 7. 



Angular Momentum 
/9=(7+5)/(2+2 7 ) 


Sound Speed 
°-=(l-7)/(2+2 7 ) 


Pressure 
A=-7/(7+l) 


Density 
1=-V(7+1) 


Disk Height 

<5=(7+3)/(2+2 7 ) 


Adiabatic Index 
7 


1.25 


-0.125 


-0.625 


-0.375 


0.875 


5/3 


1.27 


-0.115 


-0.615 


-0.385 


0.885 


8/5 


1.30 


-0.100 


-0.600 


-0.400 


0.900 


3/2 


1.33 


-0.083 


-0.583 


-0.417 


0.917 


7/5 


1.36 


-0.071 


-0.571 


-0.429 


0.929 


4/3 



Note. — These are the exponents of (r — r s ) for the various physical quantities close to the horizon; see the 
discussion in the text. 



Table 2. Model Parameters. 



Model 


7 


a 


eo 


to 


K 


K c 


r c 


Vc 


a c 


la 


r* 


V* 


a* 


1 


1.5 


0.10 


0.0 


2.60 


0.007173 


0.005222 


6.132 


0.2254 


0.2058 


2.763 


2.001 


44.6812 


0.5633 


2 


1.5 


0.30 


0.0 


1.76 


0.014750 


0.007733 


10.63 


0.1855 


0.1694 


2.258 


2.001 


44.6843 


0.6507 




1.5 


0.03 


0.0 


3.21 


0.001858 


0.001634 


4.900 


0.2084 


0.1902 


3.256 


2.001 


44.6802 


0.4299 



